Econometría II: SARIMA

Departamento de Economía

Carlos A. Yanes G.

2023-09-13

Paquetes con que se trabaja la sesión

Los paquetes que vamos a utilizar en la sesión de hoy son:

Ejecución

No olvide que debe desde luego debe tener presente en su documento de R Markdown los siguientes paquetes:

library(pacman)
p_load(readxl, TSstudio, tidyverse, stats, urca, forecast, ggfortify, ggplot2, tseries, fpp2))

Preambulo

SARIMA

  • Nos acercamos a modelos un poco mas superiores que los anteriormente vistos.

  • Llega la hora de trabajar con modelos estacionales, estos desde luego intentaran mostrar la referencia estacional en las series de tiempo.

  • Se presenta la opción de Pronostico Rápido. Se denomina auto.arima

Datos

Cartera comercial de bancos

Producto de industria Manufacturera Perú (op)

Proceso 1

Nuevamente Cartera comercial

  • La vez pasada supimos que no era una serie Ruido Blanco

  • Eso nos afecta el Pronostico de esa serie. Tal vez podamos mirar otro modelo. Puede ser estacional, como de otro orden. Eso vamos hacer hoy!!.

Cargamos datos

Code 
bd <- read_excel("cartera.xlsx")
bd <- bd |> select(Cartera)
cartera <- ts(bd, frequency=12, start=c(2016,1))
Code 
bd2 <- read_excel("Pind01.xlsx")
bd2 <- bd2 |> select(pind)
pind<- ts(bd2, frequency=12, start=c(1992,1))

Code 
cartera %>% diff(lag=12) %>% diff() %>%
  ggtsdisplay()

Code 
cartera %>% diff() %>%
  ggtsdisplay()

  • Este modelo ya diferenciado estacionalmente y -resumido (vamos paso adelante)- nos dice que el ACF cae lentamente y tenemos dos picos uno de AR y uno de MA en el lag 12 de la función de PACF. Esta parte nos indica en conformidad un ARIMA (1,1,1) para ambos correlogramas. Sin embargo los candidatos pueden ser los siguientes:

  • Los modelos que pueden ser candidatos son:

    • SARIMA \((1,1,0)(0,1,1)_{12}\)
    • SARIMA \((1,1,0)(1,1,1)_{12}\)
    • SARIMA \((1,1,1)(1,1,1)_{12}\) (opcional)

Test de Raices Unitarias

Code 
dif_cartera_12<- diff(cartera, lag = 12)
dcartera12<-diff(dif_cartera_12, 1)
dfuller<-ur.df(dcartera12, lags=0, type='trend')
summary(dfuller)

############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression trend 


Call:
lm(formula = z.diff ~ z.lag.1 + 1 + tt)

Residuals:
     Min       1Q   Median       3Q      Max 
-14554.1  -1809.3    355.4   1975.6  13263.8 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 233.54338  937.93121   0.249 0.804081    
z.lag.1      -0.39722    0.09807  -4.051 0.000129 ***
tt           -1.39875   21.84248  -0.064 0.949120    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3993 on 71 degrees of freedom
Multiple R-squared:  0.1898,    Adjusted R-squared:  0.167 
F-statistic: 8.317 on 2 and 71 DF,  p-value: 0.0005685


Value of test-statistic is: -4.0506 5.5579 8.3172 

Critical values for test statistics: 
      1pct  5pct 10pct
tau3 -4.04 -3.45 -3.15
phi2  6.50  4.88  4.16
phi3  8.73  6.49  5.47
Code 
testkp<-ur.kpss(dcartera12, type = "tau", lags = "short")
summary(testkp)

####################### 
# KPSS Unit Root Test # 
####################### 

Test is of type: tau with 3 lags. 

Value of test-statistic is: 0.083 

Critical value for a significance level of: 
                10pct  5pct 2.5pct  1pct
critical values 0.119 0.146  0.176 0.216

Code 
sarima_1<-Arima(cartera, order=c(1,1,0), seasonal=c(0,1,1))
checkresiduals(sarima_1)


    Ljung-Box test

data:  Residuals from ARIMA(1,1,0)(0,1,1)[12]
Q* = 25.588, df = 16, p-value = 0.06011

Model df: 2.   Total lags used: 18
Code 
sarima_2<-Arima(cartera, order=c(1,1,0), seasonal=c(1,1,1))
checkresiduals(sarima_2)


    Ljung-Box test

data:  Residuals from ARIMA(1,1,0)(1,1,1)[12]
Q* = 26.011, df = 15, p-value = 0.03791

Model df: 3.   Total lags used: 18
Code 
sarima_3<-Arima(cartera, order=c(1,1,1), seasonal=c(1,1,1))
checkresiduals(sarima_3)


    Ljung-Box test

data:  Residuals from ARIMA(1,1,1)(1,1,1)[12]
Q* = 18.235, df = 14, p-value = 0.1963

Model df: 4.   Total lags used: 18

  • El mejor modelo por lo pronto es el (opcional) o sarima de orden (1,1,1) en ambas partes.

  • Tenemos que igual revisar cada uno en test de ruido blanco.

Test de Ruido para cada modelo

Code 
checkresiduals(sarima_1, plot=FALSE)

    Ljung-Box test

data:  Residuals from ARIMA(1,1,0)(0,1,1)[12]
Q* = 25.588, df = 16, p-value = 0.06011

Model df: 2.   Total lags used: 18

    Ljung-Box test

data:  Residuals from ARIMA(1,1,0)(1,1,1)[12]
Q* = 26.011, df = 15, p-value = 0.03791

Model df: 3.   Total lags used: 18

    Ljung-Box test

data:  Residuals from ARIMA(1,1,1)(1,1,1)[12]
Q* = 18.235, df = 14, p-value = 0.1963

Model df: 4.   Total lags used: 18

Selección final

Usamos el criterio AKAIKE

Code 
aicc <- c(
  Arima(cartera, order=c(1,1,0), seasonal=c(0,1,1))$aicc,
  Arima(cartera, order=c(1,1,0), seasonal=c(1,1,1))$aicc,
  Arima(cartera, order=c(1,1,0), seasonal=c(1,1,0))$aicc,
  Arima(cartera, order=c(1,1,1), seasonal=c(1,1,1))$aicc
  )
  • Primer modelo SARIMA (1,1,0)(0,1,1): 1431.52.
  • Segundo modelo SARIMA (1,1,0)(1,1,1): 1433.65.
  • Tercer modelo SARIMA (1,1,0)(1,1,0): 1440.77.
  • Cuarto modelo SARIMA (1,1,1)(1,1,1): 1430.42.

Tenemos modelo 🟢

Resultado

Code 
autoplot(forecast(sarima_3, h=12))

Code 
summary(sarima_3)
Series: cartera 
ARIMA(1,1,1)(1,1,1)[12] 

Coefficients:
         ar1      ma1    sar1     sma1
      0.8453  -0.3732  0.0197  -1.0000
s.e.  0.0886   0.1373  0.1337   0.1825

sigma^2 estimated as 7503236:  log likelihood=-709.78
AIC=1429.55   AICc=1430.42   BIC=1441.14

Training set error measures:
                   ME     RMSE      MAE        MPE      MAPE       MASE
Training set 204.0135 2460.436 1628.825 0.04095422 0.3311958 0.04303489
                    ACF1
Training set -0.03445785
Code 
pb<-forecast(sarima_3, h=12)
pb
         Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
May 2023       641353.7 637602.8 645104.5 635617.3 647090.1
Jun 2023       643806.9 637132.7 650481.2 633599.5 654014.3
Jul 2023       644479.8 634796.7 654163.0 629670.8 659288.9
Aug 2023       645718.6 632979.6 658457.6 626236.0 665201.3
Sep 2023       648374.5 632577.0 664172.0 624214.3 672534.8
Oct 2023       650729.5 631903.9 669555.1 621938.2 679520.7
Nov 2023       654776.7 632976.6 676576.8 621436.3 688117.1
Dec 2023       656594.6 631889.4 681299.8 618811.2 694378.0
Jan 2024       653680.6 626150.5 681210.7 611576.9 695784.3
Feb 2024       658183.8 627925.0 688442.6 611907.0 704460.6
Mar 2024       662116.4 629207.2 695025.6 611786.2 712446.6
Apr 2024       665346.1 629866.1 700826.1 611084.2 719608.1

Auto arima

Miremos como trabaja la IA

Code 
auto.arima(cartera)
Series: cartera 
ARIMA(0,2,1)(0,0,2)[12] 

Coefficients:
          ma1    sma1    sma2
      -0.5303  0.2654  0.2823
s.e.   0.1013  0.1043  0.1207

sigma^2 estimated as 10308437:  log likelihood=-816.3
AIC=1640.59   AICc=1641.09   BIC=1650.41

Sin pasos

Code 
auto.arima(cartera, 
  stepwise=FALSE, approximation=FALSE)
Series: cartera 
ARIMA(1,2,1)(2,0,0)[12] 

Coefficients:
          ar1      ma1    sar1    sar2
      -0.0390  -0.4663  0.2833  0.2200
s.e.   0.0058   0.0121  0.0029  0.0013

sigma^2 estimated as 10191224:  log likelihood=-813.87
AIC=1637.74   AICc=1638.49   BIC=1650.01

Gracias por su atención!!

Universidad del Norte